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Abstract 

The clustering of ultrahigh energy (> 10 20 eV) cosmic rays 
(UHECR) suggests that they might be emitted by compact 
sources. Statistical analysis (Dubovsky et al., 2000) esti- 
mated the source density. We extend their analysis to give 
also the confidence intervals (CI) for the source density us- 
ing a.) no assumptions on the relationship between clustered 
and unclustered events; b.) nontrivial distributions for the 
source luminosities and energies; c.) the energy dependence 
of the propagation. We also determine the probability that a 
proton created at a distance r with energy E arrives at earth 
above a threshold E c . Using this function one can deter- 
mine the observed spectrum just by one numerical integra- 
tion for any injection spectrum. The observed 14 UHECR 
events above 10 20 eV with one doublet gives for the source 
densities 180^111° ' 10 ~ 3 M P C ~ 3 ( on the 68% confidence 
level). 

1 INTRODUCTION 

The interaction of protons with photons of the cosmic mi- 
crowave background predicts a sharp drop in the cosmic ray 
flux above the GZK cutoff around 5 • 10 19 eV (Greisen, 1966; 
Zatsepin and Kuzmin, 1966). The available data shows no 
such drop. About 20 events above 10 20 eV were observed 
by a number of experiments such as AGASA (Takeda et al., 
1998), Fly's Eye (Bird et al., 1993), Haverah Park (Lawrence 
et al., 1991), Yakutsk (Efimov et al., 1991) and HiRes (Kieda 
et al., 1999). Future experiments, particularly Pierre Auger 
(Boratav, 1996; Guerard, 1999; Bertou et al. 2000), will have 
a much higher statistics. Since above the GZK energy the 
attenuation length of particles is a few tens of megaparsecs 
(Yoshida and Teshima, 1993 ; Aharonian and Cronin, 1994; 
Protheroe and Johnson, 1996; Bhattacharjee and Sigl, 2000; 
Achterberg et al., 1999; T. Stanev et al., 2000) if an ultra- 
high energy cosmic ray (UHECR) is observed on earth it is 
usually assumed that it is produced in our vicinity. 

At these high energies the galactic and extragalactic mag- 
netic fields do not affect the orbit of the cosmic rays, thus 
they should point back to their origin within a few degrees. 
In contrast to the low energy cosmic rays one can use UHE- 
CRs for point-source search astronomy. Though there are 
some peculiar clustered events, which we discuss in detail, 



the overall distribution of UHECRs on the sky is practically 
isotropic. This observation is rather surprising since in prin- 
ciple only a few astrophysical sites are capable to accelerate 
such particles. 

There are several ways to look for the source inhomogene- 
ity from the energy spectrum and spatial directions of UHE- 
CRs. One possibility is to assume that the source density of 
UHECRs is proportional to the galaxy densities (Waxman et 
al., 1997; Giller et al., 1980; Hill and Schramm, 1985); or 
one can analyze the clustering of the unknown sources by 
some correlation length (Bahcall and Waxman 2000). 

Clearly, the arrival directions of the UHECRs measured 
by experiments show some peculiar clustering: some events 
are grouped within ~ 3°, the typical angular resolution of 
an experiment. Above 4 • 10 19 eV 92 cosmic ray events 
were detected, including 7 doublets and 2 triplets. Above 
10 20 eV one doublet out of 14 events were found (Uchihori 
et al., 2000). The chance probability of such a clustering 
from uniform distribution is rather small (Uchihori et al., 
2000; Hayashida et al., 1996; Tinyakov and Tkachev 2001a, 
Tinyakov and Tkachev 2001b,). 

The clustered features of the events initiated an interest- 
ing statistical analysis assuming compact UHECR sources 
(Dubovsky et al., 2000). The authors found a large number, 
~ 400 for the number of sources inside a GZK sphere of 
25 Mpc. We extend their analysis to give also the CIs for the 
source density using a.) no assumptions on the relationship 
between clustered and unclustered events; b.) nontrivial dis- 
tributions for the source luminosities and energies; c.) the 
energy dependence of the propagation. We also determine 
the probability that a proton created at a distance r with en- 
ergy E arrives at earth above a threshold E c . 

As we show the most probable value for the source density 
is really large; however, the statistical significance of this re- 
sult is rather weak. At present the small number of UHECR 
events allows a 95% CI for the source density which spreads 
over four orders of magnitude. Since future experiments, par- 
ticularly Pierre Auger, will have a much higher significance 
on clustering (the expected number of events of 10 20 eV and 
above is 60 per year, we present give results also for larger 
number of UHECRs above 1 20 e V. 

In order to avoid the assumptions of (Dubovsky et al., 
2000) a combined analytical and Monte-Carlo technique will 
be presented adopting the conventional picture of protons as 
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the ultrahigh energy cosmic rays. Our analytical approach 
of Section ^ gives the event clustering probabilities for any 
space, luminosity and energy distribution by using a single 
additional function P(r, E; E c ), the probability that a proton 
created at a distance r with energy E arrives at earth above 
the threshold energy E c (Bahcall and Waxman, 2000). With 
our Monte-Carlo technique of Section ^ we determine the 
probability function P(r, E; E c ) for a wide range of parame- 
ters. Our results for the present and future UHECR statistics 
are presented in Section ^[ We summarize in Section ||. 

2 ANALYTICAL APPROACH 

The number of UHECRs emitted by a source of A lumi- 
nosity during a period T follows the Poisson distribution. 
However, not all emitted UHECRs will be detected. They 
might loose their energy during propagation or can simply 
go to the wrong direction. For UHECRs the energy loss 
is dominated by the pion production in interaction with the 
cosmic microwave background radiation. In a recent anal- 
ysis (Bahcall and Waxman, 2000) the probability function 
P(r, E, E c ) was presented for three specific threshold ener- 
gies. This function gives the probability that a proton created 
at a given distance from earth (r) with some energy (E) is 
detected at earth above some energy threshold (E c ). 

The features of the Poisson distribution enforce us to take 
into account the fact that the sky is not isotropically observed. 

The probability of detecting k events from a source at dis- 
tance r with energy E can be obtained by simply including 
the factor P(r, E, E c )Aq/ (Anr 2 ) in the Poisson distribution: 

exp [-P(r, E, E c )i]j/r 2 ] 
p k (x.,E,j) = L - i- X 



is: 
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[P(r,E,E c )r)j/r 2 ] k . 
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where we introduced j = XT A/ (4tt) and Ar\j (Airr 2 ) is the 
probability that an emitted UHECR points to a detector of 
area A. The factor rj represents the visibility of the source, 
which was determined by spherical astronomy. We denote 
the space, energy and luminosity distributions of the sources 
by p(x), c(E) and h(j), respectively. The probability of de- 
tecting k events above the threshold E c from a single source 
randomly positioned within a sphere of radius R is 



P k = / dV p(x) / dE c(E) / dj h(j) x 
Js R j e c Jo 

exp [-P(r, E, E c )r,j/r 2 } r ^ & ^ 
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Denote the total number of sources within the sphere of 
sufficiently large radius (e.g. several times the GZK radius) 
by N and the number of sources that gave k detected events 
by Nk- Clearly, N = 2Vj and the total number of de- 
tected events is N e = iiVj. The probability that for N 
sources the number of different detected multiplets are N k 
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For a given set of unclustered and clustered events (iVi and 
N 2 ,N 3 ,...) inverting the P(N, {N k }) distribution function 
gives the most probable value for the number of sources and 
also the CI for it. 

Note, that P k and then P(N, {Nk}) are easily determined 
by a well behaving four-dimensional numerical integration 
for any c(E), h(j) and p(r) distribution functions. In order 
to illustrate the uncertainties and sensitivities of the results 
we used a few different choices for these distribution func- 
tions. 

For c(E) we studied three possibilities. The most straight- 
forward choice is the extrapolation of the 'conventional high 
energy component' oc E~ 2 . Another possibility is to use a 
stronger fall-off of the spectrum at energies just below the 
GZK cutoff, e.g. oc E~ 3 . The third possibility is to as- 
sume that UHECRs are some decay products of metastable 
superheavy particles (Berezinsky et al., 1997; Kuzmin and 
Rubakov, 1998; Birkel and Sarkar, 1998; Sarkar 2000; Fodor 
and Katz 2001b). According to (Birkel and Sarkar, 1998) 
the superheavy particles decay into quarks and gluons which 
initiate multi-hadron cascades through gluon bremstrahlung. 

In the recent analysis (Dubovsky et al., 2000) the authors 
have shown that for a fixed set of multiplets the minimal den- 
sity of sources can be obtained by assuming a delta-function 
distribution for h(j). We studied both this limiting lumi- 
nosity: h(j) = 5(j — and a more realistic one with 
Schechter's luminosity function (Schechter 1976), which can 
be given as: h(j)dj = h ■ (j/j*) -1 ' 25 exp(-j / 'j*)d(j /]*). 

The space distribution of sources can be given based on 
some particular survey of the distribution of nearby galax- 
ies or on a correlation length r$ characterizing the clustering 
features of sources. For simplicity the present analysis deals 
with a homogeneous distribution of sources. 

3 MONTE-CARLO STUDY OF THE 
PROPAGATION 

In our Monte-Carlo approach we determined the propaga- 
tion of UHECR protons on an event by event basis. The in- 
elasticity of Bethe-Heitler pair production is small (m 10 -3 ), 
thus we used a continuous energy loss approximation for this 
process. The inelasticity of pion-photoproduction is larger 
(« 0.2 — 0.5) in the energy range of interest, thus there are 
only a few tens of such interactions during the propagation. 
Due to the Poisson statistics and the spread of the inelastic- 
ity, we will see a spread in the energy spectrum even if the 
injected spectrum is mono-energetic. 

In our simulation protons are propagated in small steps 
(10 kpc), and after each step the energy losses due to pair pro- 
duction, pion production and the adiabatic expansion are cal- 
culated. During the simulation we keep track of the current 
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energy of the proton and its total displacement. For the pro- 
ton interaction lengths and inelasticities we used the values 
of (Bhattacharjee and Sigl, 2000; Achterberg et al., 1999). 
The deflection due to magnetic field is not taken into account 
(for a recent Monte-Carlo on it see eg. Stanev et al., 2000). 
To cover a broad energy range we used the parametrization 

P(r, E, E c ) = exp [-a • (r/1 Mpc) 6 ] . (4) 

Fig. [ljshows the functions a(E/E c ) andb(E/ E c ) forarange 
of three orders of magnitude and for five different threshold 
energies. Just using the functions of a(E/E c ) and b(E/E c ), 
thus a parametrization of P(r, E, E c ) one can obtain the ob- 
served energy spectrum for any injection spectrum without 
additional Monte-Carlo simulation. 

4 RESULTS 

In order to determine the CIs for the source densities we used 
the frequentist method (Groom et al., 2000). We wish to set 
limits on S, the source density. Using our Monte-Carlo based 
P(r, E, E c ) functions and our analytical technique we deter- 
mined p(Ni,N 2 , N 3 , S;j*), which gives the probability 
of observing Ni singlet, N 2 doublet, A3 triplet etc. events if 
the true value of the density is S and the central value of lu- 
minosity is j*. For a given set of {Ni, i — 1,2,...} the above 
probability distribution as a function of S and determines 
the 68% and 95% confidence level regions in the S—j* plane. 
For different choices of c(E) and h(j) see Table |for the cal- 
culated densities. Our results for the Dirac-delta luminosity 
distribution are in agreement with the result of (Dubovsky et 
al., 2000) within the error bars. Neverthless, there is a very 
important message. The CIs are are so large that on the 95% 
confidence level two orders of magnitude smaller densities 
than suggested as a lower bound by (Dubovsky et al., 2000) 
are also possible. 

It is of particular interest to study higher statistics too, and 
determine CIs for these cases. We performed a detailed anal- 
ysis on this question (Fodor and Katz 2001a). An interesting 
feature of the results is that "doubling" the present statistics 
with the same clustering features (in the case studied by the 
table this means one new doublet out of 10 new events) re- 
duces the CIs by an order of magnitude. The study of even 
higher statistics leads to the conclusion that experiments in 
the near future with approximately 200 UHECR events can 
tell at least the order of magnitude of the source density. 

5 SUMMARY 

We presented a technique in order to statistically analyze the 
clustering features of UHECR. The technique can be applied 
for any model of UHECR assuming small deflection. The 
key role of the analysis is played by the Pk functions defined 
by eqn. (^), which is the probability of detecting k events 
above the threshold from a single source. Using a combina- 
torial expression of eqn. (^) the probability distribution for 
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Figure 1: The functions a(E/E c ) -left panel- and 
b(E/E c ) -right panel- for the probability distribution 
function P(r, E, E c ) using the parametrization exp[— a ■ 
(r/1 Mpc) fc ] for five different threshold energies (5 • 10 19 eV, 
10 20 eV, 2 ■ 10 20 5 • 10 20 eV and 10 21 eV). 
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Table 1 : The most probable values for the source densities 
and their error bars given by the 68% and 95% confidence 
level regions (the latter in parenthesis). The numbers are in 
units of 10~ 3 Mpc~ 3 The three possible energy spectrums 
are given by a distribution proportional to E~ 2 , E~ 3 , or by 
the decay of a 10 12 GeV particle (denoted by "decay"). The 
luminosity distribution can be proportional to a Dirac-delta 
or to Schechter's luminosity function (denoted by "SLF"). 



any set of multiplets can be given as a function of the source 
density. 

We discussed several types of energy and luminosity dis- 
tributions for the sources and gave the most probable source 
densities with their CIs for present and future experiments. 

The probability P(r, E, E c ) that a proton created at a dis- 
tance r with energy E arrives above the threshold E c (Bah- 
call and Waxman, 2000) is determined and parametrized for 
a wide range of threshold energies. This result can be used 
to obtain the observed energy spectrum of the UHECR for 
arbitrary injection spectrum. 

In ref. (Dubovsky et al., 2000) the authors analyzed the 
statistical features of clustering of UHECR, which provided 
constraints on astrophysical models of UHECR if the num- 
ber of clusters is small, by giving a bound from below. In 
our paper we have shown that there is some constraint, but 
it is far from being tight. At present statistics the 95% CIs 
usually span 4 orders of magnitude. Two orders of magni- 
tude smaller numbers than the prediction of (Dubovsky et 
al., 2000) (their eqn. (13) suggests for the density of sources 
~ 6 • 10~ 3 Mpc~ 3 ) can also be obtained. Adding 10 new 
events with an additional doublet the CI can be reduced to 3 
orders of magnitude and the increase of the UHECR events 
to 200 can tell at least the order of magnitude of the source 
density. 

More details of the present analysis can be found in (Fodor 
and Katz, 2001a). Note, that a similar study based on the Z- 
burst scenario (Fargion et al., 1999; Weiler, 1999) can be car- 
ried out which gives the mass of the heaviest neutrino (Fodor 
etal.,2001). 
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